Three-body recombination of ultracold Bose gases using the truncated Wigner method 
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We apply the truncated Wigner method to the process of three-body recombination in ultracold Bose gases. 
We find that within the validity regime of the Wigner truncation for two-body scattering, three-body recombi- 
nation can be treated using a set of coupled stochastic differential equations that include diffusion terms, and 
can be simulated using known numerical methods. As an example we investigate the behaviour of a simple 
homogeneous Bose gas. 
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I. INTRODUCTION 

The dominant loss process affecting ultracold gaseous al- 
kali metal systems is inelastic three-body recombination 0, 
1^ El a process characterised by collisional events in- 
volving three atoms leading to the creation of a single two- 
atom molecule (a dimer). The binding energy released by the 
molecule formation is retained by the particles as kinetic en- 
ergy. Typically this results in the loss of all three atoms from 
the system as the molecule is not trapped by any applied ex- 
ternal potential and the energy of the free atom is high enough 
to overcome any confinement barrier. Indeed it is this pro- 
cess that limits the lifetime of experimentally produced alkali 
metal Bose-Einstein condensates, due to the large increase in 
density once the temperature is lowered past the critical point 

In a previous paper | 6] we presented a comprehensive treat- 
ment of the truncated Wigner approach for ultracold Bose 
gases including elastic two-body interactions. In this paper 
we extend that treatment to include three-body recombination 
events, which modifies the ensemble differential equations 
describing the evolution of a single realisation of the field. 
These modified differential equations are explicitly stochas- 
tic, including dynamic noise sources arising from the action of 
three-body recombination on the virtual particle background 
field. 

To provide a demonstration of our extended formalism we 
examine the evolution of a simple homogeneous system, start- 
ing from a zero-temperature state where the particle popula- 
tion is initially confined to a single (condensate) mode. 



A. Three-body recombination in ultracold gases 

Assuming that three-body recombination is the only parti- 
cle loss mechanism affecting the system, it can be shown that 
the rate of change of total particle number is 101 



dN [t] 
dt 



3K3 / dx (x, t) n (x, t)^ 



(1) 



where n (x, t) is the total number density of atoms and K3 
is the three-body recombination event rate constant. We have 
assumed that all the particles involved in the recombination 
process are lost from the system, hence the prefactor of 3 in 



Eq. which describes the number of particles lost from the 
system. The third-order normalised equiposition correlation 
function g''^^ (x, t) measures the statistics of the field, being 
unity for a fully coherent system, i.e. a zero-temperature con- 
densate, and 3! = 6 for a purely thermal system. The factor 
of 6 increase in the loss rate of thermal over coherent systems 
for similar densities has been observed experimentally |j2j] . 



II. TRUNCATED WIGNER TREATMENT 

A. The restricted field 

As in our previous work 0], we describe the many -body 
system of identical bosons using the Schrodinger picture 
bosonic field operator 



*(x)=^^,(x), 



(2) 



Here the mode operators dj annihilate a single boson from the 
jth mode, and obey the commutation relations 



0, 



a,;, a 



(3) 



while the coordinate space functions il^j (x) form an infinite 
orthonormal basis set where 



2m 



+ C/ext (x) 



i^j (x) = hwji^j (x) , 



(4) 



where C/cxt (x) is the applied external potential. 

We now divide mode space into two subspaces, a low- 
energy space (L) consisting of all those modes whose eigenen- 
ergies are less than the boundary energy ecut, and a high- 
energy space (H) that includes all remaining modes. For this 
work our interest lies with the dynamics of the low-energy 
subspace. Using these subspaces we define the field operators 



^s(x) = ^V'j(x)aj, 

j<£H 



(5a) 
(5b) 
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Of most importance for this paper is the low-energy restricted 
basis field operator '^■p (x), which can be obtained from the 
total field operator 



^(x)^ J2 V',(x), 

jeL.H 



using the projector 



(6) 



(7) 



as (x) = -P * (x) 



The restricted basis field operator 



obeys the commutation relations 



(x) , (x') 

#p(x),^';(x')j = 5p(x,x'), 

where the restricted delta function is defined by 

Sv (x,x')^EV^l(^')^'^W- 



(8) 



(9) 



The conjugate projector Q can be obtained using the comple- 
mentarity relation V + Q = 1. 



B. Master equation 



Our previous paper Ba] assumed that only two atoms par- 
ticipate in any single scattering event. In this way the particle 
interactions are described using a simple s-wave contact po- 
tential as an approximation to the full two-body T-matrix. Ob- 
viously such a description does not include three-body scatter- 
ing events. Full theoretical treatments of three-body scatter- 
ing including all possible coUisional channels are extremely 
complicated, and we do not attempt such an approach here. 
Instead we adopt a quantum-optical approach, starting from 
a phenomenologically appropriate Hamiltonian including in- 
elastic three-body recombination events, to which we apply 
the truncated Wigner method. 

We assume that within the dilute limit the charac- 
teristic range of the three-body recombination potential 
^^TBR (xi, X2, X3) is much smaller than the average interpar- 
ticle spacing. Thus, following thematically the approach for 
pairwise scattering, we replace this scattering potential by an 
effective zero-range three-body scattering T-operator, whose 
interaction strength is essentially a free parameter that will be 
chosen to satisfy experimentally observed loss rates. Within 
this approach then, in order to include the effects of three- 
body recombination the Schrodinger picture effective Hamil- 
tonian is modified to include the term 0] 



J 



H, 



(TBR) 



off 



'tTBR 



rfx 



*yx)St(x)(vl/p(x)) +(*^(x)) vI/g(x)S(x) 



(10) 



where the molecule field operator S (x) annihilates a dimer 
from the field and ktbr is a measure of the energy associated 
with the three-body process. We have assumed in formulating 
Eq. ilQ\ that the binding energy associated with the molecule 
formation is large enough that the unpaired atom generated 
by a recombination event is created within the high-energy 
subspace H, rather than the low-energy (system) subspace L, 



and is described by the high-energy field operator (x). 

Jack f?'] has considered this partial Hamiltonian, and has 
shown that by eliminating both the molecular and high-energy 
atomic fields from the evolution using a standard interaction 
picture approach for initially uncoupled fields |<S |. one obtains 
the master equation term 



dpjt) 
dt 



(TBR) 



dx 



2(*p(x)j p(t)(*;(x) 
p{t) (^t,(x))'(^p(x)) 



for the low-energy atomic subspace (system) density opera- 
tor p{t). The quantity 7 governs the rate of recombination 
events, and its relationship to /V3 we consider later. In arriv- 
ing at this master equation term it has been assumed that the 
output products of the recombination events immediately exit 
the coordinate space region containing the system, such that 
they play no further role in the evolution. To describe the full 
master equation for the system density operator one combines 



vI/T^Cx) *p(x) p{t) 



(11) 



Eq. (II 1> with the von Neumann equation calculated in 



C. Functional Wigner function correspondences 

The master equation term given by Eq. ( II U can be used 
to calculate the evolution of the corresponding multimode 
Wigner function ( { a j , a* } , t) using appropriate operator 



3 



correspondences j^l- However, rather than using the mode 
operator correspondences that were used in 0], here we per- 
form this step using functional operator correspondences. 

We define, similar to the restricted basis field operator 
^-p (x), the restricted basis wavefunctions 



VI/;, (x) ^ Y.^*{^)a] 



and the related functional derivatives 
6 



d 



(x) da, 



8^1p (x) 



da* 



(12a) 
(12b) 

(13a) 
(13b) 



Using these definitions together with the Wigner function 
mode operator correpondences f^, we find that the actions of 
the restricted basis field operator on the system density op- 
erator p (t) can be expressed as actions on the corresponding 



Wigner function using 



p(t)*p(x) ^ (*p(x)- 



1 s 


2 6"$^ 


fx) 


1 5 






(x) 


1 S 






(x) 


1 S 
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(x) 



W{t) (14a) 
VK(i) (14b) 
W{t) (14c) 
(t) .(14d) 



Such functional Wigner function operator correspondences 
have been previously used by Steel et al. ifioll . 



D. Wigner function evolution 

Applying the functional Wigner function operator corre- 
spondences, Eq. (I14> . to the master equation term describ- 
ing three-body recombination, Eq. il 1\ . we obtain, after some 
manipulation, the Wigner function evolution term 



= — ax 

dt 



( ^ 



■vl'-. 



3|vl'i,|'-9|vl'i,|'5p(x,x) + -5p(x,x)' 



9 - 18 \^vf 6v (x, x) + (x, x)' 



Am''- 



9 

+ - 



5^ 



4 V (5*^(5*^ 



53 



+ 



J5 



16 V 5^%,5'^*^ 



*pr-(5p(x,x)) +- 



16 (S*3^(5^'^3 



(15) 



Eq. (I15> is a rather complex equation of motion, including 
derivative terms up to sixth order. However, we can only 
write differential equations describing the evolution of a sin- 
gle ensemble member for Wigner function evolutions contain- 
ing derivative terms up to second order. Thus to proceed we 
must truncate the higher order terms in Eq. (11 5> , a process that 
is also required for the pairwise scattering 



1. Wigner truncation 

To justify the truncation of the higher order terms in the 
Wigner function evolution, we follow a similar method to that 
given in j6l. Let us assume that at some time r the Wigner 
function of the (low-energy) system has the factorisable form 



({a„a*}, r) ^n^cxp 



I , a," — a 



30 I 



(16) 



Here a, 



is the expectation value amplitude of the 



jth mode, and Tj is proportional to the inverse width of the 



Wigner function for that mode. This type of function de- 
scribes both coherent (where Vj = 2) and thermally dis- 
tributed modes, but does not describe number states or other 
more exotic states. The factorisability of this Wigner function 
indicates that number fluctuations between disparate modes 
are uncorrected. 

Evaluating the Wigner function evolution given by Eq. O 
using the Wigner function given by Eq. returns a rather 
complicated expression, which we give in full in the appendix, 
Eq. (lAU . Essentially we find that for increasing order in 
the leading order term in |vl/p| decreases. In per- 
forming the Wigner truncation for the two-body scattering 
in i6j], we required that in the coordinate space regions of 
high real particle density that n (x) ^ 5-p (x, x). Given that 
I '5 73 1 ^ n and that all remaining terms scale as 5-p, and as- 
suming that there is significant real particle density in the re- 
gions where three-body recombination is important compared 
to the local density of modefunctions, only the first few terms 
in Y^-p I are important. By keeping only those terms of 4th and 
5th order in \ we find that, within the same validity regime 
for the two-body elastic scattering, the Wigner function equa- 
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tion of motion can be accurately described by 



dt 



dx 



r 



J2 



(17) 



While it is possible to directly convert this equation of mo- 
tion into a set of coupled differential equations, the functional 
nature of the derivative operators can obscure some of the 
details. Instead we choose to perform the conversion using 



J 



an explicit mode representation. Using our definitions of the 
functional derivatives, Eq. il3\ . we find that the Wigner func- 
tion evolution due to three-body recombination can be ex- 
pressed as 



^(TBR) 

dt 



E 



d 7 



W 



(18) 



r 



E. Stochastic differential equations 

It is important to remember when converting Eq. (I18> to 
its equivalent differential equations that we have two sets of 
independent variables, {aj} and {a^}- Thus while the drift 
terms are straightforward, and we find by using the relations 
given in fllll for the Ito calculus that 



■ J dxi'jl'^SJvf-^v. (19) 



where Aj<- = A* as required, the diffusion terms are not so 
easily obtained. 



To obtain the terms in the stochastic differential equations 
corresponding to the diffusion terms in Eq. ( I18> . we first find 
it necessary to rewrite the coefficient of that diffusive part as 



dx / dx'V^;i*p|' E i'kii^'.n^'vfi'l 

keL.H 

^ / dxiJ*\^j,Uk f dx' ii^'.r i^'T^f 



(20) 
(21) 

(22) 



keL,H 



r 



where it is important to note that the summation over the index 
k runs over the complete mode space L ® H. Including basis 
modes that are not part of the system subspace into the for- 
malism in this way may appear to be cause for concern, as the 
master equation term from which we are working, Eq. (Ill> 
contains no reference to these high-energy states. However, 
we note that we are free to choose any set of ensemble dif- 
ferential equations that can be shown to be mathematically 
equivalent to the Fokker-Planck equation lITlll . such that our 



inclusion of the high-energy modes in Eq. 
mathematically accurate. 



22b is certainly 



It can be shown, either by rewriting Eq. ( I18> in terms of ex- 
plicitly real quantities (including the mode amplitude quadra- 
tures) and directly using the conversion relations given in 
or by working backwards using complex Ito calculus, that the 
ensemble differential equations corresponding to the diffusive 
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part of the Wigner function evolution are given by 



which can be straightforwardly shown to obey 



k£L.H 

for all those modes j E L. Here the complex Wiener pro- 
cesses dWk (t) obey the relations 

(dWkit)) = (24a) 
{dWkit)dWi(t)) = (24b) 
{dWk (t) dW; (t)) = Skjdt. (24c) 



{dW{x,t)) = 
{dW {Tc,t)dW (x',t)) = 
{dW{yi,t)dW*{x',t)) = (5(x- 



x') dt. 



(26a) 
(26b) 
(26c) 



In fact, given the local nature of the recombination process 
in coordinate space, a more useful form of the total Wiener 
process is given by 

dWi^,t)= i^k{^)dWkit), (25) 

keLM 



Inserting the spatial Wiener process dW as appropriate into 
the diffusive mode evolution given by Eq. (I23> . using the drift 
mode evolutions of Eq. ( I19> and including the evolution in the 
absence of three-body recombination given in |6], gives the 
total evolution of the low-energy system mode amplitudes as 



daj = —iujjaj dt + / dxt/j* < — 



^■pdt+ ^^\-^-pfdW (x,t)l. 



(27) 



Using the definition of the system wavefunction given by Eq. il2\ we find that the corresponding evolution of the coordinate 
space field is 



2m 



^rdt + V 



^>T,dt+\j^\^>-p\^dW{^,t)Y (28) 



r 



where we have recognised the low-energy projector V, 
Eq. 0. 



1. Rate of population change 
The total (real) particle population of the field is defined as 

= > (N,)=y ■ (d,a) . (29) 



Using the correspondence of moments of the Wigner function 
to symmetrically ordered products of quantum operators 
we find that N can be calculated using Wigner function aver- 
ages as 



TV : 



E 



1 



(30) 



w 



Using this result, we find the rate of change of particle number 
for a single trajectory to be given by Ito's formula 



dN 
It 



lim — (da*ai + a*dai + da*dai) . (31) 



The terms | da* daj } are included here because, unlike ordi- 
nary deterministic calculus, the presence of the Wiener pro- 
cesses in Eq. (127 > give these terms a non-zero value in the 
limit dt ^ 0. 



Taking expectation values, using the properties of the 
spatially-dependent Wiener process, Eq. (I26> . we find the en- 
semble averaged rate of normalisation change to be 



dN 
'dt 



-7 



dx 



w 



w 



(32) 

where we have written 6-p for S-p (x, x) for compactness, as 
we also do below. It may appear from Eq. (I32t that our trun- 
cated Wigner treatment of three-body recombination intro- 
duces a small correction to the rate of particle loss, apparently 
creating particles (in the average). However, a clearer under- 
standing can be obtained by expressing the moments of the 
Wigner function as physically significant quantities. 

Using the properties of the Wigner function moments we 
find, for example, that 



2(*|,*p)<5p + -<5^(33) 



—S-p, 
2 ' ' 



(34) 



where we have again suppressed the spatial dependences. Re- 
placing the Wigner function moments in Eq. ( I32> in this way 
returns 

'dN' 

g'- 'W -t- '71 O-p T -Zll-'Jp 



dt 



-1 



dx 



(35) 



Thus, rather than reducing the rate of particle loss, the trun- 
cated Wigner treatment leads to an increased rate of particle 



6 



loss. However, given that we have required that n ^ 5-p 
to perform the Wigner truncation, this correction should be 
small. Note that Eq. (I35> also shows that particle loss only oc- 
curs in those regions where there is real particle density, such 
that those coordinate space regions solely occupied by virtual 
particles will exhibit zero particle loss. 

Comparing Eq. ( I35> to Eq. Q shows that 7 = iK^. Thus 
while A'3 is the rate constant for three-body recombination 
events, 7 is the number loss rate constant. 

It is worthwhile discussing a possible point of confusion 
when using these classical field methods. As the field for a 
single trajectory is represented by a single wavefunction, it 
could be considered that the field is therefore uniformly co- 
herent at all points. In such a case those behaviours that de- 
pend upon the statistics of the field, such as three-body recom- 
bination, would be improperly treated. However, this view 
is incorrect, as such statistics only obtain physical meaning 
when considering ensembles of trajectories. As an example, 
while direct inclusion of the statistics is relevant when consid- 
ering three-body recombination using the mean (either spa- 
tial or temporal) particle density, in those regions where the 
system exhibits thermal statistics, the trajectory wavefunction 
will exhibit density fluctuations both spatially and temporally. 
Thus the (spatial and temporal) mean of \'$-pf will be larger 
than the mean density cubed, leading to the increased rate of 
loss observed experimentally. Indeed, given that a single tra- 
jectory is entirely analogous to a single experimental run, the 
fact that here the wavefunction contains the full behaviour of 
the field is, if not obvious, at least eminently reasonable. This 
result also provides for the factor of two increase in nonlinear 
interaction strength between the condensate and the thermal 
particles due to the exchange energy lll2il . 

F. Plane wave basis 

While our formalism is applicable to any orthonormal 
single-particle basis {ipj (x)}, the most useful set of mode- 
functions for many situations, including the simple system we 
consider in this paper, is the plane-wave modes. For this basis 
the modes are eigenstates of the curvature operator only, such 
that Uc^t = and hujj = ti^kj /2m where kj = |kj |, and 

^^{^) = ±=e^^'-\ (36) 
V V 

Within a periodically bounded volume of extent V = y< 
Ly X the orthonormal plane-wave modes are arranged in 
momentum space such that 

k,=^k. + ^k.^k., (37) 

^x 

where mj, rij and pj are integers. Using this plane-wave ba- 
sis, the energy cutoff that defines the low-energy mode sub- 
space becomes a spherical cutoff in momentum space, with 
the boundary defined by ftfccut = V2me^. 

III. NUMERICAL SIMULATIONS 

To demonstrate our truncated Wigner treatment of three- 
body recombination we have numerically simulated a 



(relatively) simple zero-temperature homogeneous gas of 
mp) — |1, —1) ■^'^Na atoms. We describe the system us- 
ing a set of plane-wave modes, with the initial real particle 
population confined to the ground (kj = 0) mode. 

Determination of the three-body recombination event rate 
constant for various alkali metals has been performed both 
theoretically QUI [3 and experimentally ISBUlIll- For 
this paper we take as a best estimate of the relevant the 
value measured at MIT fl4ll for a fully condensed gas 

A3 = 3.7 X 10"^^ cm^s"^ (38) 

In that work it was reported that optical confinement was 
used to produce a rather large particle density of ^'^Na of 
3 X 10^^ cm^"^ at the centre of the trap. Thus, given that 
the rate of particle loss scales as and the correction due to 
the dynamic noise sources as n?, such a large particle den- 
sity should provide information on a parameter regime where 
three-body recombination is significant. 

We use a simulation volume of = (4.7/Ltm)'^, such that to 
achieve an initial (uniform) density of n = 3 x 10^^ cm^"^ we 
use iVo (i = 0) = 3.06 x 10^. The mode spacing (in velocity 
space) along each of the cartesian directions, Eq. (13 7> . is de- 
termined by the volume to be 3.7 mms^^. To characterise the 
strength of the interactions, we use a = 2.75 nm. 

We have performed simulations using two distinct bound- 
aries to the low-energy subspace; v^nt = 44.3 mms^^, for 
which the number of modes M — 7.2 x 10'^; and Wcut = 
59.1 mms~^, for which M = 1.7 x 10^. In both cases the 
number of modes is significantly less than the number of real 
particles, and we therefore expect that both cutoffs will return 
valid results. 



A. Initial states 

For our zero-temperature homogeneous system, the appro- 
priate initial state for a single trajectory is described by 

ao (0) = [^0 + iBo] aj#o (0) = ^ [A, + iB^] . 

(39) 

Here Aj and Bj are Gaussian random variables of zero mean 
and unit variance, such that 

{A,) = {B,) = {A.,B,) = 
{A,Aj) = {B,B,)^1. (40) 

This initial state, Eq. ( I39> satisfies the assumed Wigner func- 
tion used to justify the Wigner truncation, Eq. il6\ . with 

Tj = 2 for all j, ao^ = \/lVo and aj^Og = 0. 

B. Evolution algorithm 

The dynamic noise term present in the Eq. ( I27> means that 
we cannot directly apply the deterministic projected RK4IP 
algorithm, which was used to obtain the results of 01 ■ Rather 
we must employ an algorithm that explicitly allows for such 
time-dependent random processes. 

The simplest such method is the Euler algorithm j3, in 
which the drift terms are calculated at the start of each time 
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step and the continuous Wiener processes dW are replaced by 
a single discrete Wiener process AW. However, for any rea- 
sonable accuracy the time step for an Euler algorithm must be 
very small, thus requiring very long calculation times. Mil- 
stein and Tretyakov liQ] have considered various more so- 
phisticated algorithms for propagating stochastic differential 
equations with dynamic random processes. In particular, they 
have shown that for situations where the influence of the dy- 
namic noise on the system is very much less than the deter- 
ministic evolution (the small noise limit), one can accurately 
describe the total evolution using a relatively simple modifi- 
cation to the fourth-order Runge-Kutta (RK4) algorithm. Es- 
sentially, in this method one calculates the deterministic evo- 
lution using the RK4 algorithm, while the dynamic noise is 
calculated using an Euler type derivative calculation based on 
the state of the system at the start of each time step. This re- 
sult therefore allows us to use a slightly modified version of 
the projected RK4IP algorithm to propagate Eq. ( I27l l. 

Importantly, any numerical propagation method requires a 
discrete coordinate space, such that the relations given for the 
spatial Wiener process, Eq. j26l do not apply. Rather, we use 
noise sources that obey 



particle number 



{dW,, (t) dW, (t)) 

{dw^it)dw:it)) 






1 

aV 



(41a) 
(41b) 

(41c) 



where dW^ {t) is the time-dependent Wiener process at the 
/ith point on the coordinate space simulation grid and AV is 
the volume space increment about that grid point. The algo- 
rithm advances the field in time by the increment At with each 
application, and we use At = 250 ns for all our simulations. 

The total number of real particles within the low-energy 
subspace is given by Eq. OOt . where the subtraction of 1/2 
can be understood as removing the virtual particles introduced 
into the initial state of the field, Eq. ( I39> . For systems with a 
large number of modes M, an excellent estimate of the total 
particle number can be made using 



N{t) 



El" 



(42) 



In Fig. n we plot the estimated total (real) particle numbers 
calculated using Eq. (I42> for single trajectories of the sys- 
tem described above using cutoffs of Vcut = 44.3 mms~^ 
and 59.1 mms^^. From these curves we observe that, over 
the larger time-scale, the total particle populations of the sys- 
tems decrease, apparently mono tonic ally, with the Vcut = 
59.1 mms^^ trajectory showing greater particle loss. On the 
smaller time-scale however, as shown by the inset, the total 
particle populations fluctuate rapidly, on a scale of roughly 
1-10 particles per time-step. 

To provide a comparison with our truncated Wigner results, 
consider a simple model. For a homogeneous system, and 
assuming that the third-order correlation function is both spa- 
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(3) 



Nit) 



N{Q) 



6K3g(3)N(0f 





(43) 
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tially and temporally invariant, such that g'^' (x, t) 
Eq. Q can be integrated to return the time-dependent total sions. 



FIG. 1: Total particle numbers for the sample trajectories of the 
system described in the text. Main plot shows the populations for 
(higher and lower solid lines) Vcut ~ {44.3, 59.1} mms^^, together 
with the results for the simple model (dashed line). The inset shows 
the population for the trajectory with Vcut = 44.3 mms~^ over a 
subset of the total time. 



For our system, using g^^' = 1, the model shows a smaller 
rate of loss than that observed for either of our trajectories, as 
shown by the dashed line in Fig.^ This result is predicted by 
Eq. ( I35> . as is the difference in the two trajectory populations. 



C. Results 

IV. CONCLUSION 

The truncated Wigner description of ultracold Bose gases 
has many significant advantages over more traditional ap- 
proaches, such as the Gross-Pitaevskii equation, and the ex- 
tension outlined in this paper allows for the inclusion of three- 
body recombination processes. We have shown that within the 
validity regime of the Wigner truncation for two-body scatter- 
ing, three-body recombination can be described using stochas- 
tic differential equations describing the evolution of a single 
trajectory, which can be solved using numerical techniques. 
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APPENDIX A: MATHEMATICAL DETAILS 



Using the particular Wigner function given by Eq. (I16> in the full three-body recombination Wigner function equation of 
motion, Eq. il5i . and evaluating all d/daj, d/da* operators returns 

'jTi^(TBR) r ( / 



(3rd) 



(4th) 



(6th) 



J6i 



J6i 



[(e; - *P + (Cp - enJ] 6 - I' - 6 E Y 1^^- 1' 



J6i 



6iep-^7.oi'-3EYi^^i'r^ 



J6i 



(5th) - [(ef - *p + *^ (Cp - cpo)] icp - I' - E y 



'jeL 



3 ICp - ?P„r - 6 - enJ' E Y l^^l' + M E Y l^^-l' 



S-p 



|iep-e7.or-3ie7.-eni'EYi^^i' 



J6i 



-3iep-SPoiMEYi^^i') -(EyI'/^^i' 



> W. 



r 



(Ala) 



(Alb) 



(Ale) 



(Aid) 



(Ale) 



(Alf) 



Here the terms arising from the separate derivative orders are derivative terms are given on the first and second lines). For 
kept separated, as indicated on the left-hand side (first order compactness we have written S-p = Sp (x, x). 
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